For the total Lagrangian approach, the discrete equations are formulated with respect to the reference configuration. For the updated Lagrangian approach, the discrete equations are formulated in the current configuration , which is assumed to be the new reference configuration .
The independent variables in the total Lagrangian approach are
X
{\displaystyle X}
and
t
{\displaystyle t}
. In the updated Lagrangian they are
x
{\displaystyle x}
and
t
{\displaystyle t}
which are with respect to the new reference configuration.
The dependent variable in the total Lagrangian approach is the displacement
u
(
X
,
t
)
{\displaystyle u(X,t)}
. In the updated Lagrangian approach, the dependent variables are the Cauchy stress
σ
(
x
,
t
)
{\displaystyle \sigma (x,t)}
and the velocity
v
(
x
,
t
)
{\displaystyle v(x,t)}
.
The stress measure is the Cauchy stress :
σ
(
x
,
t
)
=
T
A
.
{\displaystyle {\sigma (x,t)={\cfrac {T}{A}}~.}}
The strain measure is the rate of deformation (also called velocity strain):
D
(
x
,
t
)
=
∂
v
∂
x
=
v
,
x
.
{\displaystyle {D(x,t)={\frac {\partial v}{\partial x}}=v_{,x}~.}}
Note that the derivative is a spatial derivative. It can be shown that
∫
0
t
D
(
x
,
τ
)
d
τ
=
ln
(
F
(
x
,
t
)
)
{\displaystyle \int _{0}^{t}D(x,\tau )~d\tau =\ln(F(x,t))}
where
F
{\displaystyle F}
is the deformation gradient. So the integral of the rate of deformation in 1-D is similar to the logarithmic strain (also called natural strain).
The Updated Lagrangian governing equations are:
ρ
J
=
ρ
0
.
{\displaystyle {\rho ~J=\rho _{0}~.}}
For the rod,
ρ
A
F
=
ρ
0
A
0
.
{\displaystyle \rho ~A~F=\rho _{0}~A_{0}~.}
These are the same as those for the total Lagrangian approach.
∂
∂
x
(
A
σ
)
+
ρ
A
b
=
ρ
A
0
D
v
D
t
.
{\displaystyle {{\frac {\partial }{\partial x}}(A~\sigma )+\rho ~A~b=\rho ~A_{0}~{\cfrac {Dv}{Dt}}~.}}
In this case
A
{\displaystyle A}
may not be constant at along the length and
further simplification cannot be done. In short form,
(
A
σ
)
,
x
+
ρ
A
b
=
ρ
A
v
˙
.
{\displaystyle (A\sigma )_{,x}+\rho ~A~b=\rho ~A~{\dot {v}}~.}
If we ignore heat flux and heat sources
ρ
∂
w
int
∂
t
=
σ
D
.
{\displaystyle {\rho {\frac {\partial w_{\text{int}}}{\partial t}}=\sigma ~D~.}}
If we include heat flux and heat sources
ρ
∂
w
int
∂
t
=
σ
D
−
∂
q
∂
x
+
ρ
s
.
{\displaystyle {\rho {\frac {\partial w_{\text{int}}}{\partial t}}=\sigma ~D-{\frac {\partial q}{\partial x}}+\rho ~s~.}}
In short form:
ρ
w
˙
int
=
σ
D
−
q
,
x
+
ρ
s
.
{\displaystyle \rho {\dot {w}}_{\text{int}}=\sigma ~D-q_{,x}+\rho ~s~.}
For a linear elastic material:
σ
(
x
,
t
)
=
E
σ
D
∫
0
t
D
(
x
,
τ
)
d
τ
=
E
σ
D
ln
(
F
(
x
,
t
)
)
.
{\displaystyle {\sigma (x,t)=E^{\sigma D}\int _{0}^{t}D(x,\tau )~d\tau =E^{\sigma D}\ln(F(x,t))~.}}
The superscript
σ
D
{\displaystyle \sigma D}
refers to the fact that this function
relates
σ
{\displaystyle \sigma }
and
D
{\displaystyle D}
.
For small strains:
E
σ
D
=
Youngs modulus
.
{\displaystyle E^{\sigma D}=~~{\text{Youngs modulus}}~.}
For a linear elastic material:
σ
˙
(
x
,
t
)
=
E
σ
D
D
(
x
,
t
)
.
{\displaystyle {{\dot {\sigma }}(x,t)=E^{\sigma D}~D(x,t)~.}}
For the updated Lagrangian approach, initial conditions are needed
for the stress and the velocity. The initial displacement is assumed
to be zero.
The initial conditions are:
σ
(
x
,
0
)
=
σ
0
(
x
)
for
x
∈
[
0
,
L
]
v
(
x
,
0
)
=
v
0
(
x
)
for
x
∈
[
0
,
L
]
u
(
x
,
0
)
=
0
for
x
∈
[
0
,
L
]
{\displaystyle {\begin{aligned}\sigma (x,0)&=\sigma _{0}(x)&~{\text{for}}~&x\in [0,L]\\v(x,0)&=v_{0}(x)&~{\text{for}}~&x\in [0,L]\\u(x,0)&=0&~{\text{for}}~&x\in [0,L]\end{aligned}}}
The essential boundary conditions are
v
(
x
,
t
)
=
v
¯
(
x
,
t
)
for
x
∈
Γ
v
.
{\displaystyle {v(x,t)={\bar {v}}(x,t)\qquad ~{\text{for}}~x\in \Gamma _{v}~.}}
The traction boundary conditions are
n
σ
(
x
,
t
)
=
t
x
(
x
,
t
)
for
x
∈
Γ
t
.
{\displaystyle {n~\sigma (x,t)=t_{x}(x,t)\qquad ~{\text{for}}~x\in \Gamma _{t}~.}}
The unit normal to the boundary in the current configuration is
n
{\displaystyle n}
. The tractions in the current configuration are related to those
in the reference configuration by
t
x
A
=
t
x
0
A
0
.
{\displaystyle {t_{x}~A=t_{x}^{0}~A_{0}~.}}
The interior continuity or jump condition is
⌊
σ
A
⌉
=
0
.
{\displaystyle {\lfloor \sigma ~A\rceil =0~.}}
The momentum equation is
(
A
σ
)
,
x
+
ρ
A
b
=
ρ
A
D
v
D
t
.
{\displaystyle (A~\sigma )_{,x}+\rho ~A~b=\rho ~A~{\cfrac {Dv}{Dt}}~.}
To get the weak form over an element, we multiply the equation by a
weighting function and integrate over the current length of the
element (from
x
a
{\displaystyle x_{a}}
to
x
b
{\displaystyle x_{b}}
).
∫
x
a
x
b
w
[
(
A
σ
)
,
x
+
ρ
A
b
−
ρ
A
D
v
D
t
]
d
x
=
0
.
{\displaystyle \int _{x_{a}}^{x_{b}}w\left[(A~\sigma )_{,x}+\rho ~A~b-\rho ~A~{\cfrac {Dv}{Dt}}\right]~dx=0~.}
Integrate the first term by parts to get
∫
x
a
x
b
w
(
A
σ
)
,
x
d
x
=
[
w
A
σ
]
x
a
x
b
−
∫
x
a
x
b
w
,
x
(
A
σ
)
d
x
.
{\displaystyle \int _{x_{a}}^{x_{b}}w(A~\sigma )_{,x}~dx=\left[w~A~\sigma \right]_{x_{a}}^{x_{b}}-\int _{x_{a}}^{x_{b}}w_{,x}~(A~\sigma )~dx~.}
Plug the above into the weak form and rearrange to get
∫
x
a
x
b
w
,
x
(
A
σ
)
d
x
+
∫
x
a
x
b
w
ρ
A
D
v
D
t
d
x
=
∫
x
a
x
b
w
ρ
A
b
d
x
+
[
w
A
σ
]
x
a
x
b
.
{\displaystyle \int _{x_{a}}^{x_{b}}w_{,x}~(A~\sigma )~dx+\int _{x_{a}}^{x_{b}}w~\rho ~A~{\cfrac {Dv}{Dt}}~dx=\int _{x_{a}}^{x_{b}}w~\rho ~A~b~dx+\left[w~A~\sigma \right]_{x_{a}}^{x_{b}}~.}
If we think of the weighting function
w
{\displaystyle w}
as a variation of
v
{\displaystyle v}
that
satisfies the kinematic admissibility conditions, we get
the principle of virtual power :
∫
x
a
x
b
(
δ
v
)
,
x
(
A
σ
)
d
x
+
∫
x
a
x
b
δ
v
ρ
A
D
v
D
t
d
x
=
∫
x
a
x
b
δ
v
ρ
A
b
d
x
+
[
δ
v
A
σ
]
x
a
x
b
.
{\displaystyle {\int _{x_{a}}^{x_{b}}(\delta v)_{,x}(A~\sigma )~dx+\int _{x_{a}}^{x_{b}}\delta v~\rho ~A~{\cfrac {Dv}{Dt}}~dx=\int _{x_{a}}^{x_{b}}\delta v~\rho ~A~b~dx+\left[\delta v~A~\sigma \right]_{x_{a}}^{x_{b}}~.}}
Recall that
v
,
x
=
D
(
x
,
t
)
and
t
x
=
n
σ
.
{\displaystyle v_{,x}=D(x,t)\qquad {\text{and}}\qquad t_{x}=n~\sigma ~.}
Therefore,
δ
v
,
x
=
δ
D
(
x
,
t
)
{\displaystyle \delta v_{,x}=\delta D(x,t)}
and
[
δ
v
A
σ
]
x
a
x
b
=
[
δ
v
A
t
x
]
Γ
t
.
{\displaystyle \left[\delta v~A~\sigma \right]_{x_{a}}^{x_{b}}=\left[\delta v~A~t_{x}\right]_{\Gamma _{t}}~.}
Hence we can alternatively write the weak form as
∫
x
a
x
b
δ
D
A
σ
d
x
+
∫
x
a
x
b
δ
v
ρ
A
D
v
D
t
d
x
=
∫
x
a
x
b
δ
v
ρ
A
b
d
x
+
[
δ
v
A
t
x
]
Γ
t
.
{\displaystyle {\int _{x_{a}}^{x_{b}}\delta D~A~\sigma ~dx+\int _{x_{a}}^{x_{b}}\delta v~\rho ~A~{\cfrac {Dv}{Dt}}~dx=\int _{x_{a}}^{x_{b}}\delta v~\rho ~A~b~dx+\left[\delta v~A~t_{x}\right]_{\Gamma _{t}}~.}}
Comparing with the energy equation, we see that the first term above is
the internal virtual power. Using the substitutions
A
d
x
=
d
Ω
{\displaystyle A~dx=d\Omega }
and
D
v
/
D
T
=
v
˙
{\displaystyle Dv/DT={\dot {v}}}
, we can also write the weak form as
∫
Ω
δ
D
σ
d
Ω
+
∫
Ω
δ
v
ρ
v
˙
d
Ω
=
∫
Ω
δ
v
ρ
b
d
Ω
+
[
δ
v
A
t
x
]
Γ
t
.
{\displaystyle {\int _{\Omega }\delta D~\sigma ~d\Omega +\int _{\Omega }\delta v~\rho ~{\dot {v}}~d\Omega =\int _{\Omega }\delta v~\rho ~b~d\Omega +\left[\delta v~A~t_{x}\right]_{\Gamma _{t}}~.}}
The weak form may also be written in terms of the virtual powers as
δ
P
=
δ
P
int
−
δ
P
ext
+
δ
P
kin
=
0
{\displaystyle {\delta P=\delta P_{\text{int}}-\delta P_{\text{ext}}+\delta P_{\text{kin}}=0}}
where,
δ
P
int
=
∫
x
a
x
b
(
δ
v
)
,
X
(
A
σ
)
d
x
δ
P
ext
=
∫
x
a
x
b
δ
v
ρ
A
b
d
x
+
[
δ
v
A
t
X
]
Γ
t
δ
P
kin
=
∫
x
a
x
b
δ
v
ρ
A
v
˙
d
x
{\displaystyle {\begin{aligned}\delta P_{\text{int}}&=\int _{x_{a}}^{x_{b}}(\delta v)_{,X}~(A~\sigma )~dx\\\delta P_{\text{ext}}&=\int _{x_{a}}^{x_{b}}\delta v~\rho ~A~b~dx+\left[\delta v~A~t_{X}\right]_{\Gamma _{t}}\\\delta P_{\text{kin}}&=\int _{x_{a}}^{x_{b}}\delta v~\rho ~A~{\dot {v}}~dx\end{aligned}}}
The trial velocity field is
v
h
(
x
,
t
)
=
∑
j
=
1
n
N
j
(
x
)
v
j
(
t
)
{\displaystyle {v_{h}(x,t)=\sum _{j=1}^{n}N_{j}(x)v_{j}(t)}}
In matrix form,
v
h
(
x
,
t
)
=
N
(
x
)
v
(
t
)
where
N
=
[
N
1
(
X
)
N
2
(
X
)
N
n
(
X
)
]
and
v
=
[
v
1
(
t
)
v
2
(
t
)
⋮
v
n
(
t
)
]
.
{\displaystyle {v_{h}(x,t)=\mathbf {N} (x)~\mathbf {v} (t)}~~{\text{where}}~~\mathbf {N} ={\begin{bmatrix}N_{1}(X)&N_{2}(X)&&N_{n}(X)\end{bmatrix}}~{\text{and}}~\mathbf {v} ={\begin{bmatrix}v_{1}(t)\\v_{2}(t)\\\vdots \\v_{n}(t)\end{bmatrix}}~.}
The resulting acceleration field is the material time derivative
of the velocity
a
h
(
x
,
t
)
=
v
˙
h
(
x
,
t
)
=
∑
j
=
1
n
N
j
(
x
)
v
˙
j
(
t
)
.
{\displaystyle a_{h}(x,t)={\dot {v}}_{h}(x,t)=\sum _{j=1}^{n}N_{j}(x){\dot {v}}_{j}(t)~.}
Hence,
a
h
(
x
,
t
)
=
∑
j
=
1
n
N
j
(
x
)
a
j
(
t
)
or
a
h
(
X
,
t
)
=
N
(
x
)
a
(
t
)
where
a
=
[
a
1
(
t
)
a
2
(
t
)
⋮
a
n
(
t
)
]
.
{\displaystyle {a_{h}(x,t)=\sum _{j=1}^{n}N_{j}(x)a_{j}(t)}~{\text{or}}~{a_{h}(X,t)=\mathbf {N} (x)~\mathbf {a} (t)}~~{\text{where}}~~\mathbf {a} ={\begin{bmatrix}a_{1}(t)\\a_{2}(t)\\\vdots \\a_{n}(t)\end{bmatrix}}~.}
Note that the shape functions are functions of
X
{\displaystyle X}
and not of
x
{\displaystyle x}
. We
could transform them into function of
x
{\displaystyle x}
using the inverse mapping
x
=
φ
−
1
(
x
,
t
)
.
{\displaystyle x=\varphi ^{-1}(x,t)~.}
However, in that case the shape functions become functions of time and
the procedure becomes more complicated.
The test (weighting) function is
δ
v
h
(
x
)
=
∑
i
=
1
n
N
i
(
x
)
δ
v
i
.
or
δ
v
h
(
x
)
=
N
δ
v
where
δ
v
=
[
δ
v
1
δ
v
2
⋮
δ
v
n
]
.
{\displaystyle {\delta v_{h}(x)=\sum _{i=1}^{n}N_{i}(x)\delta v_{i}~.}~~{\text{or}}~~{\delta v_{h}(x)=\mathbf {N} ~\delta \mathbf {v} }~~{\text{where}}~~\delta \mathbf {v} ={\begin{bmatrix}\delta v_{1}\\\delta v_{2}\\\vdots \\\delta v_{n}\end{bmatrix}}~.}
The derivatives of the test functions with respect to
x
{\displaystyle x}
are
(
δ
v
h
)
,
x
=
∑
i
=
1
n
N
i
,
x
δ
v
i
.
or
(
δ
v
h
)
,
x
=
B
δ
v
where
B
=
[
N
1
,
x
N
2
,
x
N
n
,
x
]
.
{\displaystyle {(\delta v_{h})_{,x}=\sum _{i=1}^{n}N_{i,x}\delta v_{i}~.}~~{\text{or}}~~{(\delta v_{h})_{,x}=\mathbf {B} ~\delta \mathbf {v} }~~{\text{where}}~~\mathbf {B} ={\begin{bmatrix}N_{1,x}&N_{2,x}&&N_{n,x}\end{bmatrix}}~.}
It is convenient to use the isoparametric concept to compute the spatial
derivatives. Let us reexamine what this approach involves.
Figure 1 shows a two-noded 1-D element in the reference and
current configurations along with the parent element.
Figure 1. Reference and Current Configurations of a 1-D element.
The map between the Eulerian coordinates and the parent coordinates is
x
(
ξ
,
t
)
=
(
1
−
ξ
)
x
1
e
(
t
)
+
ξ
x
2
e
(
t
)
ξ
=
N
1
(
ξ
)
x
1
e
(
t
)
+
N
2
(
ξ
)
x
2
e
(
t
)
or
x
(
ξ
,
t
)
=
N
(
ξ
)
x
e
(
t
)
.
{\displaystyle x(\xi ,t)=(1-\xi )~x_{1}^{e}(t)+\xi ~x_{2}^{e}(t)~\xi =N_{1}(\xi )~x_{1}^{e}(t)+N_{2}(\xi )~x_{2}^{e}(t)\qquad {\text{or}}~~{x(\xi ,t)=\mathbf {N} (\xi )~\mathbf {x} ^{e}(t)~.}}
At
t
=
0
{\displaystyle t=0}
,
x
(
ξ
,
0
)
=
x
(
ξ
)
=
N
(
ξ
)
X
e
.
{\displaystyle x(\xi ,0)=x(\xi )=\mathbf {N} (\xi )~\mathbf {X} ^{e}~.}
Therefore the displacement in the parent coordinates is
u
(
ξ
,
t
)
=
x
(
ξ
,
t
)
−
x
(
ξ
)
=
N
(
ξ
)
[
x
e
(
t
)
−
X
e
]
or
u
(
ξ
,
t
)
=
N
(
ξ
)
u
e
(
t
)
.
{\displaystyle u(\xi ,t)=x(\xi ,t)-x(\xi )=\mathbf {N} (\xi )~\left[\mathbf {x} ^{e}(t)-\mathbf {X} ^{e}\right]\qquad {\text{or}}~~{u(\xi ,t)=\mathbf {N} (\xi )~\mathbf {u} ^{e}(t)~.}}
Similarly, the velocity and its variation in the parent coordinates are
given by
v
(
ξ
,
t
)
=
N
(
ξ
)
v
e
(
t
)
,
δ
v
(
ξ
,
t
)
=
N
(
ξ
)
δ
v
e
(
t
)
.
{\displaystyle {v(\xi ,t)=\mathbf {N} (\xi )~\mathbf {v} ^{e}(t)~,\qquad \delta v(\xi ,t)=\mathbf {N} (\xi )~\delta \mathbf {v} ^{e}(t)~.}}
The acceleration in the parent coordinates is given by
a
(
ξ
,
t
)
=
N
(
ξ
)
v
˙
e
(
t
)
.
{\displaystyle {a(\xi ,t)=\mathbf {N} (\xi )~{\dot {\mathbf {v} }}^{e}(t)~.}}
The rate of deformation is given by
D
(
x
,
t
)
=
v
,
x
(
x
,
t
)
=
N
,
x
(
X
(
x
,
t
)
)
v
e
(
t
)
.
{\displaystyle D(x,t)=v_{,x}(x,t)=\mathbf {N} _{,x}(X(x,t))~\mathbf {v} ^{e}(t)~.}
We can evaluate the derivative with respect to
x
{\displaystyle x}
by simply using the
map to the parent coordinate instead of mapping back to the reference
coordinates. Using the chain rule (for the two noded element)
N
1
,
ξ
=
N
1
,
x
x
,
ξ
and
N
2
,
ξ
=
N
2
,
x
x
,
ξ
.
{\displaystyle N_{1,\xi }=N_{1,x}~x_{,\xi }\qquad {\text{and}}\qquad N_{2,\xi }=N_{2,x}~x_{,\xi }~.}
In matrix form,
N
,
ξ
=
N
,
x
x
,
ξ
⟹
N
,
x
=
N
,
ξ
(
x
,
ξ
)
−
1
=:
B
(
ξ
)
.
{\displaystyle {\mathbf {N} _{,\xi }=\mathbf {N} _{,x}~x_{,\xi }}\qquad \implies \qquad {\mathbf {N} _{,x}=\mathbf {N} _{,\xi }~\left(x_{,\xi }\right)^{-1}=:\mathbf {B} (\xi )~.}}
The rate of deformation in parent coordinates is then given by
D
(
ξ
,
t
)
=
B
(
ξ
)
v
e
(
t
)
.
{\displaystyle {D(\xi ,t)=\mathbf {B} (\xi )~\mathbf {v} ^{e}(t)~.}}
To derive the finite element system of equations, we follow the usual
approach of substituting the trial and test functions into the
approximate weak form
∫
x
a
x
b
(
δ
v
h
)
,
x
(
A
σ
)
d
x
+
∫
x
a
x
b
δ
v
h
ρ
A
D
v
h
D
t
d
x
=
∫
x
a
x
b
δ
v
h
ρ
A
b
d
x
+
[
δ
v
h
A
σ
]
x
a
x
b
.
{\displaystyle \int _{x_{a}}^{x_{b}}(\delta v_{h})_{,x}(A~\sigma )~dx+\int _{x_{a}}^{x_{b}}\delta v_{h}~\rho ~A~{\cfrac {Dv_{h}}{Dt}}~dx=\int _{x_{a}}^{x_{b}}\delta v_{h}~\rho ~A~b~dx+\left[\delta v_{h}~A~\sigma \right]_{x_{a}}^{x_{b}}~.}
Let us proceed term by term.
The first term represents the virtual internal power
δ
P
int
=
∫
x
a
x
b
(
δ
v
h
)
,
x
(
A
σ
)
d
x
.
{\displaystyle \delta P_{\text{int}}=\int _{x_{a}}^{x_{b}}(\delta v_{h})_{,x}(A~\sigma )~dx~.}
Plugging in the derivative of the test function, we get
δ
P
int
=
∫
x
a
x
b
[
∑
i
=
1
n
N
i
,
x
δ
v
i
]
(
A
σ
)
d
x
=
∑
i
=
1
n
δ
v
i
[
∫
x
a
x
b
N
i
,
x
(
A
σ
)
d
x
]
=
∑
i
=
1
n
δ
v
i
f
i
int
.
{\displaystyle {\begin{aligned}\delta P_{\text{int}}&=\int _{x_{a}}^{x_{b}}\left[\sum _{i=1}^{n}N_{i,x}\delta v_{i}\right](A~\sigma )~dx\\&=\sum _{i=1}^{n}\delta v_{i}\left[\int _{x_{a}}^{x_{b}}N_{i,x}(A~\sigma )~dx\right]\\&=\sum _{i=1}^{n}\delta v_{i}f_{i}^{\text{int}}~.\end{aligned}}}
The matrix form of the expression for virtual internal work is
δ
P
int
=
δ
v
T
f
int
where
f
int
=
[
f
1
int
f
2
int
⋮
f
n
int
]
.
{\displaystyle {\delta P_{\text{int}}=\delta \mathbf {v} ^{T}\mathbf {f} _{\text{int}}}~~{\text{where}}~~\mathbf {f} _{\text{int}}={\begin{bmatrix}f_{1}^{\text{int}}\\f_{2}^{\text{int}}\\\vdots \\f_{n}^{\text{int}}\end{bmatrix}}~.}
The internal force is
f
i
int
=
∫
x
a
x
b
N
i
,
x
(
A
σ
)
d
x
=
∫
−
1
1
N
i
,
ξ
(
x
,
ξ
)
−
1
(
A
σ
)
x
,
ξ
d
ξ
=
∫
−
1
1
N
i
,
ξ
(
A
σ
)
d
ξ
.
{\displaystyle {f_{i}^{\text{int}}=\int _{x_{a}}^{x_{b}}N_{i,x}(A~\sigma )~dx=\int _{-1}^{1}N_{i,\xi }(x_{,\xi })^{-1}(A~\sigma )~x_{,\xi }~d\xi =\int _{-1}^{1}N_{i,\xi }(A~\sigma )~d\xi ~.}}
Note that the above simplification only occurs in 1-D.
In matrix form,
f
int
=
∫
x
a
x
b
B
T
(
A
σ
)
d
x
=
∫
−
1
1
(
N
,
ξ
)
T
(
A
σ
)
d
ξ
.
or
f
int
=
∫
Ω
B
T
σ
d
Ω
.
{\displaystyle {\mathbf {f} _{\text{int}}=\int _{x_{a}}^{x_{b}}\mathbf {B} ^{T}~(A~\sigma )~dx=\int _{-1}^{1}(\mathbf {N} _{,\xi })^{T}~(A~\sigma )~d\xi ~.}~~{\text{or}}~~{\mathbf {f} _{\text{int}}=\int _{\Omega }\mathbf {B} ^{T}~\sigma ~d\Omega ~.}}
Remark:
Note that we can write
N
,
x
=
N
,
X
d
X
d
x
⟹
N
,
x
d
x
=
N
,
X
d
X
.
{\displaystyle N_{,x}=N_{,X}~{\cfrac {dX}{dx}}\qquad \implies \qquad N_{,x}~dx=N_{,X}~dX~.}
If we use the above relation, we get
f
int
=
∫
x
a
x
b
(
N
,
x
)
T
σ
A
d
x
=
∫
X
a
X
b
(
N
,
X
)
T
σ
A
d
X
.
{\displaystyle \mathbf {f} _{\text{int}}=\int _{x_{a}}^{x_{b}}(\mathbf {N} _{,x})^{T}~\sigma ~A~dx=\int _{X_{a}}^{X_{b}}(\mathbf {N} _{,X})^{T}~\sigma ~A~dX~.}
Using the relation
σ
A
=
P
A
0
{\displaystyle \sigma ~A=P~A_{0}}
we get
f
int
=
∫
X
a
X
b
(
N
,
X
)
T
P
A
0
d
X
.
{\displaystyle {\mathbf {f} _{\text{int}}=\int _{X_{a}}^{X_{b}}(\mathbf {N} _{,X})^{T}~P~A_{0}~dX~.}}
The above is the same as the expression we had for the internal force
in the total Lagrangian formulation.
The second term represents the virtual kinetic power
δ
P
kin
=
∫
x
a
x
b
δ
v
h
ρ
A
D
v
h
D
t
d
x
.
{\displaystyle \delta P_{\text{kin}}=\int _{x_{a}}^{x_{b}}\delta v_{h}~\rho ~A~{\cfrac {Dv_{h}}{Dt}}~dx~.}
Plugging in the test function, we get
δ
P
kin
=
∫
x
a
x
b
[
∑
i
=
1
n
δ
v
i
N
i
]
ρ
A
D
v
h
D
t
d
x
=
∑
i
=
1
n
δ
v
i
[
∫
x
a
x
b
ρ
A
N
i
D
v
h
D
t
d
x
]
=
∑
i
=
1
n
δ
v
i
f
i
kin
.
{\displaystyle {\begin{aligned}\delta P_{\text{kin}}&=\int _{x_{a}}^{x_{b}}\left[\sum _{i=1}^{n}\delta v_{i}N_{i}\right]\rho ~A~{\cfrac {Dv_{h}}{Dt}}~dx\\&=\sum _{i=1}^{n}\delta v_{i}\left[\int _{x_{a}}^{x_{b}}\rho ~A~N_{i}~{\cfrac {Dv_{h}}{Dt}}~dx\right]\\&=\sum _{i=1}^{n}\delta v_{i}f_{i}^{\text{kin}}~.\end{aligned}}}
The inertial (kinetic) force is
f
i
kin
=
∫
x
a
x
b
ρ
A
N
i
D
v
h
D
t
d
x
.
{\displaystyle {f_{i}^{\text{kin}}=\int _{x_{a}}^{x_{b}}\rho ~A~N_{i}~{\cfrac {Dv_{h}}{Dt}}~dx~.}}
Now, plugging in the trial function into the expression for
the inertial force, we get
f
i
kin
=
∫
x
a
x
b
ρ
A
N
i
[
∑
j
=
1
n
D
v
j
D
t
N
j
]
d
x
=
∑
j
=
1
n
[
∫
x
a
x
b
ρ
A
N
i
N
j
d
x
]
D
v
j
D
t
=
∑
j
=
1
n
M
i
j
a
j
{\displaystyle {\begin{aligned}f_{i}^{\text{kin}}&=\int _{x_{a}}^{x_{b}}\rho ~A~N_{i}\left[\sum _{j=1}^{n}{\cfrac {Dv_{j}}{Dt}}N_{j}\right]~dx\\&=\sum _{j=1}^{n}\left[\int _{x_{a}}^{x_{b}}\rho ~A~N_{i}~N_{j}~dx\right]{\cfrac {Dv_{j}}{Dt}}\\&=\sum _{j=1}^{n}M_{ij}a_{j}\end{aligned}}}
where
M
i
j
:=
∫
x
a
x
b
ρ
A
N
i
N
j
d
x
←
Consistent Mass Matrix Coefficient.
{\displaystyle {M_{ij}:=\int _{x_{a}}^{x_{b}}\rho ~A~N_{i}~N_{j}~dx\qquad \leftarrow \quad {\text{Consistent Mass Matrix Coefficient.}}}}
and
a
j
:=
D
v
j
D
t
=
v
˙
j
←
Acceleration.
{\displaystyle a_{j}:={\cfrac {Dv_{j}}{Dt}}={\dot {v}}_{j}\qquad \leftarrow \quad {\text{Acceleration.}}}
In matrix form,
f
kin
=
M
a
where
a
=
[
a
1
a
2
⋮
a
n
]
=
[
v
˙
1
v
˙
2
⋮
v
˙
n
]
and
f
kin
=
[
f
1
kin
f
2
kin
⋮
f
n
kin
]
.
{\displaystyle {\mathbf {f} _{\text{kin}}=\mathbf {M} \mathbf {a} }~~{\text{where}}~~\mathbf {a} ={\begin{bmatrix}a_{1}\\a_{2}\\\vdots \\a_{n}\end{bmatrix}}={\begin{bmatrix}{\dot {v}}_{1}\\{\dot {v}}_{2}\\\vdots \\{\dot {v}}_{n}\end{bmatrix}}~~~{\text{and}}~~\mathbf {f} _{\text{kin}}={\begin{bmatrix}f_{1}^{\text{kin}}\\f_{2}^{\text{kin}}\\\vdots \\f_{n}^{\text{kin}}\end{bmatrix}}~.}
The consistent mass matrix in matrix notation is
M
=
∫
x
a
x
b
ρ
A
N
T
N
d
x
=
∫
Ω
ρ
N
T
N
d
Ω
.
{\displaystyle {\mathbf {M} =\int _{x_{a}}^{x_{b}}\rho ~A~\mathbf {N} ^{T}~\mathbf {N} ~dx=\int _{\Omega }\rho ~\mathbf {N} ^{T}~\mathbf {N} ~d\Omega ~.}}
Plugging the expression for the inertial force into the expression
for virtual kinetic power, we get
δ
P
kin
=
∑
i
=
1
n
δ
v
i
[
∑
j
=
1
n
M
i
j
a
j
]
=
∑
i
=
1
n
∑
j
=
1
n
δ
v
i
M
i
j
a
j
{\displaystyle {\begin{aligned}\delta P_{\text{kin}}&=\sum _{i=1}^{n}\delta v_{i}\left[\sum _{j=1}^{n}M_{ij}a_{j}\right]\\&=\sum _{i=1}^{n}\sum _{j=1}^{n}\delta v_{i}M_{ij}a_{j}\end{aligned}}}
In matrix form,
δ
P
kin
=
(
δ
v
)
T
M
a
=
(
δ
v
)
T
f
kin
.
{\displaystyle {\delta P_{\text{kin}}=(\delta \mathbf {v} )^{T}\mathbf {M} \mathbf {a} =(\delta \mathbf {v} )^{T}~\mathbf {f} _{\text{kin}}~.}}
Remark:
Note that since the integration domain is a function of time, the mass
matrix for the updated Lagrangian formulation is also a function of time.
However, from the conservation of mass, we have
ρ
0
A
0
d
X
=
ρ
A
d
x
.
{\displaystyle \rho _{0}~A_{0}~dX=\rho ~A~dx~.}
Therefore, we can write the mass matrix as
M
=
∫
x
a
x
b
ρ
A
N
T
N
d
x
=
∫
X
a
X
b
ρ
0
A
0
N
T
N
d
X
=
∫
Ω
0
ρ
0
N
T
N
d
Ω
0
.
{\displaystyle {\mathbf {M} =\int _{x_{a}}^{x_{b}}\rho ~A~\mathbf {N} ^{T}~\mathbf {N} ~dx=\int _{X_{a}}^{X_{b}}\rho _{0}~A_{0}~\mathbf {N} ^{T}~\mathbf {N} ~dX=\int _{\Omega _{0}}\rho _{0}~\mathbf {N} ^{T}~\mathbf {N} ~d\Omega _{0}~.}}
Hence the mass matrices are the same for both formulations.
The right hand side terms represent the virtual external power
δ
P
ext
=
∫
x
a
x
b
δ
v
h
ρ
A
b
d
x
+
[
δ
v
h
A
t
x
]
Γ
t
.
{\displaystyle \delta P_{\text{ext}}=\int _{x_{a}}^{x_{b}}\delta v_{h}~\rho ~A~b~dx+\left[\delta v_{h}~A~t_{x}\right]_{\Gamma _{t}}~.}
Plugging in the test function into the above expression gives
δ
P
ext
=
∫
x
a
x
b
[
∑
i
=
1
n
N
i
δ
v
i
]
ρ
A
b
d
x
+
[
(
∑
i
=
1
n
N
i
δ
v
i
)
A
t
x
]
Γ
t
=
∑
i
=
1
n
δ
v
i
[
∫
x
a
x
b
N
i
ρ
A
b
d
x
]
+
[
∑
i
=
1
n
δ
v
i
N
i
A
t
X
]
Γ
t
=
∑
i
=
1
n
δ
v
i
[
∫
x
a
x
b
N
i
ρ
A
b
d
x
+
[
N
i
A
t
x
]
Γ
t
]
=
∑
i
=
1
n
δ
v
i
f
i
ext
.
{\displaystyle {\begin{aligned}\delta P_{\text{ext}}&=\int _{x_{a}}^{x_{b}}\left[\sum _{i=1}^{n}N_{i}~\delta v_{i}\right]\rho ~A~b~dx+\left[\left(\sum _{i=1}^{n}N_{i}~\delta v_{i}\right)A~t_{x}\right]_{\Gamma _{t}}\\&=\sum _{i=1}^{n}\delta v_{i}\left[\int _{x_{a}}^{x_{b}}N_{i}~\rho ~A~b~dx\right]+\left[\sum _{i=1}^{n}\delta v_{i}~N_{i}~A~t_{X}\right]_{\Gamma _{t}}\\&=\sum _{i=1}^{n}\delta v_{i}\left[\int _{x_{a}}^{x_{b}}N_{i}~\rho ~A~b~dx+\left[N_{i}~A~t_{x}\right]_{\Gamma _{t}}\right]\\&=\sum _{i=1}^{n}\delta v_{i}~f_{i}^{\text{ext}}~.\end{aligned}}}
In matrix notation,
δ
P
ext
=
δ
v
T
f
ext
where
f
ext
=
[
f
1
ext
f
2
ext
⋮
f
n
ext
]
.
{\displaystyle {\delta P_{\text{ext}}=\delta \mathbf {v} ^{T}~\mathbf {f} _{\text{ext}}}~~{\text{where}}~~\mathbf {f} _{\text{ext}}={\begin{bmatrix}f_{1}^{\text{ext}}\\f_{2}^{\text{ext}}\\\vdots \\f_{n}^{\text{ext}}\end{bmatrix}}~.}
The external force is given by
f
i
ext
=
∫
x
a
x
b
N
i
ρ
A
b
d
x
+
[
N
i
A
t
x
]
Γ
t
.
{\displaystyle {f_{i}^{\text{ext}}=\int _{x_{a}}^{x_{b}}N_{i}~\rho ~A~b~dx+\left[N_{i}~A~t_{x}\right]_{\Gamma _{t}}~.}}
In matrix notation,
f
ext
=
∫
x
a
x
b
N
T
ρ
A
b
d
x
+
[
N
T
A
t
x
]
Γ
t
or
f
ext
=
∫
Ω
ρ
N
T
b
d
Ω
+
[
N
T
A
t
x
]
Γ
t
.
{\displaystyle {\mathbf {f} _{\text{ext}}=\int _{x_{a}}^{x_{b}}\mathbf {N} ^{T}~\rho ~A~b~dx+\left[\mathbf {N} ^{T}~A~t_{x}\right]_{\Gamma _{t}}}~~{\text{or}}~~{\mathbf {f} _{\text{ext}}=\int _{\Omega }\rho ~\mathbf {N} ^{T}~b~d\Omega +\left[\mathbf {N} ^{T}~A~t_{x}\right]_{\Gamma _{t}}~.}}
Remark:
Using the conservation of mass
ρ
0
A
0
d
x
=
ρ
A
d
x
{\displaystyle \rho _{0}~A_{0}~dx=\rho ~A~dx}
and the relations between the traction sin the reference and the current
configurations
A
0
t
x
0
=
A
t
x
{\displaystyle A_{0}~t_{x}^{0}=A~t_{x}}
we can transform the integral in the expression for the
external force to one over the reference coordinates as follows:
f
ext
=
∫
x
a
x
b
N
T
ρ
0
A
0
b
d
x
+
[
N
T
A
0
t
x
0
]
Γ
t
.
{\displaystyle \mathbf {f} _{\text{ext}}=\int _{x_{a}}^{x_{b}}\mathbf {N} ^{T}~\rho _{0}~A_{0}~b~dx+\left[\mathbf {N} ^{T}~A_{0}~t_{x}^{0}\right]_{\Gamma _{t}}~.}
The above is the same as the expression we had for the external force
in the total Lagrangian formulation.
The finite element equations in updated Lagrangian form are:
v
(
x
,
t
)
=
∑
i
N
i
(
x
)
v
i
e
(
t
)
=
N
v
e
D
(
x
,
t
)
=
∑
i
N
i
,
x
v
i
e
(
t
)
=
B
v
e
f
int
e
=
∫
Ω
e
B
T
σ
d
Ω
f
ext
e
=
∫
Ω
e
ρ
N
T
b
d
Ω
+
[
N
T
A
t
x
]
Γ
t
e
M
e
=
∫
Ω
0
e
ρ
0
N
T
N
d
Ω
0
M
u
¨
=
f
ext
−
f
int
.
{\displaystyle {\begin{aligned}v(x,t)&=\sum _{i}N_{i}(x)~v_{i}^{e}(t)=\mathbf {N} ~\mathbf {v} ^{e}\\D(x,t)&=\sum _{i}N_{i,x}~v_{i}^{e}(t)=\mathbf {B} ~\mathbf {v} ^{e}\\\mathbf {f} _{\text{int}}^{e}&=\int _{\Omega ^{e}}\mathbf {B} ^{T}~\sigma ~d\Omega \\\mathbf {f} _{\text{ext}}^{e}&=\int _{\Omega ^{e}}\rho ~\mathbf {N} ^{T}~b~d\Omega +\left[\mathbf {N} ^{T}~A~t_{x}\right]_{\Gamma _{t}^{e}}\\\mathbf {M} ^{e}&=\int _{\Omega _{0}^{e}}~\rho _{0}~\mathbf {N} ^{T}~\mathbf {N} ~d\Omega _{0}\\\mathbf {M} ~{\ddot {\mathbf {u} }}&=\mathbf {f} _{\text{ext}}-\mathbf {f} _{\text{int}}~.\end{aligned}}}